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Abstract 

A numerical program is presented which facihtates a computation pertaining to 
the fuh set of one-gluon loop diagrams (including ghost loop contributions), with M 
attached external gluon lines in all possible ways. The feasibility of such a task rests on 
a suitably defined master formula, which is expressed in terms of a set of Grassmann 
and a set of Feynman parameters. The program carries out the Grassmann integration 
and performs the Lorentz trace on the involved functions, expressing the result as 
a compact sum of parametric integrals. The computation is based on tracing the 
structure of the final result, thus avoiding all intermediate unnecessary calculations 
and directly writing the output. Similar terms entering the final result are grouped 
together. The running time of the program demonstrates its effectiveness, especially 
for large M. 
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PROGRAM SUMMARY 

Manuscript title: A numerical algorithm for efficient computations of one-gluon loop 
Feynman diagrams in QCD for a large number of external gluons 
Authors: A. S. Kapoyannis, A. I. Karanikas and C. N. Ktorides 
Program title: DIL0G2 
Programming language: FORTRAN 90 

Computer (s) for which the program has been designed: Personal Computer 

Operating system(s) for which the program has been designed: Windows 98, XP, LINUX 

Number of processors used: one 

Keywords: one-gluon loop, Feynman diagram, QCD diagram 
PACS: 02.70.Rw, 12.38.Bx 

CPC Library Classification: 11.5 Quantum Chromodynamics, Lattice Gauge Theory 

External routines/libraries used: none 

CPC Program Library subprograms used: none 

Nature of problem: The computation of one gluon/ghost loop diagrams in QCD with 
many external gluon lines is a time consuming task, practically beyond reasonable reach of 
analytic procedures. We apply recently proposed master formulas towards the computation 
of such diagrams with an arbitrary number (M) of external gluon lines, achieving a final 
result which reduces the problem to one involving integrals over the standard set, for given 
M, of Feynman parameters. 

Solution method: The structure of the master expressions is analysed from a numerical 
computation point of view. Using the properties of Grassmann variables we identify all the 
different forms of terms that appear in the final result. Each form is called "structure". We 
calculate theoretically the number of terms belonging to every "structure" . We carry out the 
calculation organising the whole procedure into separate calculations of the terms belonging 
to every "structure" . Terms which do not contribute to the final result are thereby avoided. 
The final result, extending to large values of M, is also presented with terms belonging to 
the same "structure" grouped together. 

Restrictions: M is coded as a 2-digit integer. Overflow in the dimension of used array is 
expected to appear for M > 20 in a processor that uses 4-bytes integers or for M > 34 in a 
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processor with 8-bytes integers. 

Running time: Depends on M, see enclosed figures. 
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LONG WRITE-UP 
1. Introduction 

Quantum Chromodynamics (QCD) enjoys universal acceptance as the fundamental the- 
ory for the strong interaction. As a quantum field theoretical system, QCD has been exten- 
sively applied to situations in which its perturbative content provides a dependable compu- 
tational tool. It is, in fact, within the framework of this perturbative content that QCD has 
successfully confronted the quantitative description of the multitude of scattering processes, 
which probe strong interaction dynamics at high energies. Admittedly, the study of the 
non-perturbative domain of the theory offers intriguing and, most certainly, fundamental 
challenges. Nevertheless, the immediate need to confront recent measurements coming from 
the HERA and Tevatron particle accelerators as well as the expected ones, in the near future, 
from the LHC accelerator continues to put perturbative QCD (pQCD) at the forefront of 
theoretical activity. 

Given the non-abelian structure of QCD, the (by far) most demanding component of the 
theory, with respect to perturbative calculations, is its gluonic, as opposed to its quark, sec- 
tor^. In particular, perturbative computations involving Feynman diagrams with gluon/ghost 
loops become, to say the least, quite monstrous. During the last decade or so various meth- 
ods, aiming to expedite Feynman diagram computations in QCD, have been proposed whose 
basic feature is that they rely in a first, rather than the usual second, quantization approach 
to the formulation of the theory. Corresponding attempts have employed either strings [1-3] , 
or world-line paths [4-11] as their underlying basic agents. Within the framework of the 
latter case, two of the present authors [10,11], have been involved in work which led to the 
formulation of a set of master expressions, that condense the multitude of all Feynman di- 
agrams entering a given configuration. These expressions are determined by the number of 
loops and the number of external gluon propagators attached to them. To be more precise, 
the derived expressions go up to two loop configurations, nevertheless the "logic" of the con- 
struction can be extended to loops of higher order. Suffice it to say, at this point, that the 
analytical confrontation of a two loop situation with four "external" gluon lines constitutes 
a challenging enough problem [12]. 

^For that matter, this is more so the case for the non-perturbative domain of the theory. 
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The basic feature of the master expressions arrived at in [10,11] is that they are furnished 
in terms of a set of Grassman and a set of Feynman variables. Once integrations over these 
two sets of variables arc performed one obtains the full result, i.e. the one which, for the 
given configuration, contains the contribution of all Feynman diagrams at once. It is obvious, 
even before laying an eye on these master formulas, that in order to put them into practical 
use, their confrontation calls for the employment of suitable computational methods. It is 
the aim of this paper to present a program, which will be applied to the one gluon loop 
case for M external gluonic lines (M fairly large). The main part of our program deals with 
the confrontation of multi-Grassmann variable integrals and arrives at expressions which 
involve the appropriate set of Feynman parameters only. This program could hopefully find 
applicability to other situations, where multi-Grassmann variable integrals also make their 
entrance. 

Our paper is organized as follows. In the following section we present the battery of 
formulas, which are associated with the master expression corresponding to one gluon/ghost 
loop with M external gluon attachments in all possible ways. We shall consider cases up 
to fairly large values of M, which exhibit divergent terms^, in addition to finite ones. In 
section 3 we present the ideas on which this program is based along with the formulas which 
count the terms appearing in our final result. In section 4 we describe the structure of the 
program, while section 5 presents our results, accompanied by tables and figures. Finally, 
our concluding remarks are made in section 6. 

2. The one loop master formula 

Consider a configuration consisting of one gluon/ghost loop onto which M external gluon 
lines, with corresponding momenta pi ■ ■ -pM are attached (see Figure 1). According to [10], 
the master expression, which summarizes the total contribution from all Feynman diagrams 
pertaining to this configuration is given by 

^As expected, the aforementioned master expressions implicate the absence of divergent terms for M > 4, 
cf. [10]. 
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where g is the couphng constant of the theory, the tQ , i = 1, . . . , M are the SU {S) color group 
generators (in the adjoint representation) with Trc the trace over the color group, the Ui 
are Feynman parameters, the function 9 is specified by 



9{um, ...,ui)^ 9{um - um-i) ■ ■ ■ 9{u2 - ui) 



and 



f(^)K,...,«m;T) 



n JdUd^n (Tri$W-2) 



.n=M ' 



X 



X exp 



■ (2) 



In the above relation the ^'s are Grassmann variables, the are polarization vectors for 
the external gluons, $'^1 is the so-called spin factor entering the world-line description of 
QCD (see below), with Tri denoting trace with respect to Lorenz generator representation 
indices and the G{un,Um) arc free propagators for the particle modes entering the worldline 
path integral description of QCD, in the context of its first quantized version (see [10]). They 
obey the equation (s) 



- dndmG{Un,Um) = d'^G{Un,Um) = G{Un,Um) = 2[S{Un,Um) - 1] , 



(3) 



with boundary condition 



dnG{Un, Urn) = G{Un, Um) = sign{Un - Um) - 2{Un - Um) = -G{Um, Un) ■ (4) 

The explicit expression for the spin factor in terms of the set of parameters entering our 
expressions is (the J^^, are the Lorentz generators, in the vector representation) 



= Pexp 



M 
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where 

A point of note is the following: In the above expressions a specific time ordering has 
been chosen according to which index n + 1 comes immediately after index n, with ^m+i = 
— 0. 

The saturation of indices p, a in Jp^ is performed instantly, since 

(V2)(J,.)^.2(£>,. - = -2{elPn,. - elPn,,) , (7) 

miJp,),.{4/T){e;^h: - e^+h;) = -{A/T){e;-^h: - e^+h;) (8) 
Thus, eq. (5) may be rewritten as 

M M n2-l 

= S,. + (-1) E '^.^ W + (-1)' E E <l>,^in2)Mn^) + .... (9) 

n=l n2=l ni=l 

3. Theoretical considerations surrounding the computation 

The object of computation is the quantity F^^\ as given by eq. (2). The main effort 
amounts to carrying out the Grassmann integrations entering this expression. In between 
there intervenes the task of performing trace operations over Lorentz indices associated with 
the loop(s). The tracing involves strings of the (j)^uin), cf. eq. (6). The Lorentz trace can 
be postponed until after the Grassmann integration. This trace accounts for setting the 
first and the last Lorentz index in a series of products of the objects (pnvip) equal. Since 
the second index of a (t>ij,u{n) factor must be saturated with first index of the factor that 
follows, the Grassmann integrations will be carried out first if the specific order by which 
the sequence of the (p^yinys are placed in each product is kept undisturbed. Accordingly, 
the (pfj^uin) are considered non-commutative objects during the Grassmann integration. By 
activating this rule the Lorentz indices can be dropped until after the Grassmann integration 
has been completed. 

To get a concrete handle on the situation, we introduce the quantities 

^„ = -2(£X--^>n,M), (10) 
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and 
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Using the fact that for commutative objects exp(C + D) = exp(C) exp{D), the Grassmann 
calculation assumes the form 

1 
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A direct way to proceed with this calculation is to suitably code the functions and the 
Grassmann variables and then separately calculate the objects exp CnCnC'„j 



exp 



Y^mitn ^nCnCm^mDnm ^-Ild | P exp 



— 2|. Each of these 



Sn(CnCn^Ti + ^n+l^n+l^n^nBn) 

objects may contain sums of products with at most M Grassmann variables Then the 
multiplication between them can be carried out. Since the following Grassmann properties 
are valid 

J d^n^n 1) J d^n^n 1) J d^n J d^ji 0, ^nCrra CmCn) CnCm CmCn; (1^) 

the output of the integration becomes obvious. Only those terms survive which contain 
exactly M products of ^n^n with all the indices different from each other. This calculation 
is straightforward and we have written out the corresponding FORTRAN code for carrying 
it out [13]. The routine is effective (we present in Fig. 3 results up to M = 10), but the time 
consumed to reach the result rises exponentially with larger values of M. The basic obstacle 
to the effectiveness of the routine is that the plethora of terms involved in the intermediate 
computations does not survive in the final result. A second source of trouble is the consumed 
amount of time for rearranging terms at each stage of the computation before identical ones 
are gathered and summed. 

In this paper the aforementioned issues are confronted by following a different compu- 
tational procedure, which is based on the following objectives: (a) at each stage of the 



calculation only the surviving terms are computed and (b) similar terms are grouped to- 
gether in the final result. The above goals become more concrete by making the following 
observations concerning the form of the terms entering our final expressions. These terms 
will, in general, be products of the quantities A„, C„ and Dnm accompanied by Grass- 
mann variables the origin of each one of which becomes evident from eq. (14). It is also 
apparent that each of An or C„ is accompanied by one pair of Grassmann variables, whereas 
two such pairs accompany Bn and Dnm- Let A^i, A^2, -^3 and A^4 be the number of the 
functions C„ and Dnmi respectively appearing in a given term of the final result. 

According to the Grassmann properties, the terms entering this result must satisfy 

Ni + 2N2 + N^ + 2Ni = M , < A^. < M, TVs + ^"4 7^ 1, A", g N, i = 1, 2, 3, 4 , (16) 

with the condition A^3 + A4 7^ 1 coming from the fact that the Lorentz trace of a single 
function An or Bn is zero. 

The above equation is very important because it helps one to determine all the possible 
forms of the final terms before actually carrying out the calculation. Every group of the 
integers A, represents a particular form of the final terms, each of which will be referred to 
as "structure". Finding all the possible groups of Aj suffices to determine all the possible 
structures. In this way the calculation is broken into a number of subcalculations in each of 
which one collects all the terms that correspond to a specific structure. A term belonging 
to a particular structure cannot be similar to a term belonging to another structure, since 
such terms will contain different numbers of functions Am Bn, Cn and Dnm, so goal (b) is 
satisfied. 

Let us turn now to goal (a). Since in the final result there must exist exactly M pairs 
of Grassmann variables, all different from each other, it is evident that every part of one 
exponential entering eq. (14) has to be multiplied by those components of the rest of the 
exponentials which satisfy cq. (16). We shall try now to trace how many terms survive in a 
specific structure characterised by the integers Aj. 

Firstly, as far as the Pexp function entering (14) is concerned, one may ask how can 
one obtain exactly A^i functions An and A2 functions S„ placed in an ordered product 
according to a specific time order. It is better to count the terms of the substructure 
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A^i , N2 first because the fact that the associated product is ordered imposes more restrictive 
conditions. For every group of Grassmann pairs taken from the total M pairs which enter 
the substructure N1.N2. the remaining Grassmann pairs will always form a substructure 
N3, N4 (according to (16)). The converse, however, is not true: If a substructure A^3, A^4 is 
formed first, then the remaining Grassmann pairs cannot always form a substructure with 
given N2. The reason is that the quantities are associated with the Grassmann factor 
^n+iCn+i^n^n, counccting adjaccut Grassmann pairs. Consider, for example, the case where 
A''2 = 1 and the time order is of the form 1, 2, 3, • ■ •. If, now, we are left with the Grassmann 
indices 1,3, after considering the previous terms, then it is obvious that such a case will not 
lead to surviving terms. 

/ iVi + A^2 ^ 

Generally speaking there are ways one can have products of Ni func- 

V ) 

tions and N2 Bn functions, if one is not concerned with the index n they carry. This result 
counts the ways the Ni An functions (considered indistinguishable) can occupy + N2 
total vacant places. The rest are then filled with the functions -B„. It can be checked that 
for everyone of these products of functions the ways the indices n can be distributed among 
the functions of the same kind is the same. So let us consider a product where the first N2 
functions are of the kind Bn (contributing 2N2 Grassmann pairs) and the last A^i functions 
are of the kind An (contributing A^^i Grassmann pairs). If the first 2N2 Grassmann pairs 
(according to the given time order) are attributed to the A^2 functions then there are 



^ M-2N2 



ways this 



M — 2N2 Grassmann pairs left for the A^i functions An and there are 

can be achieved. Let us now suppose that the first 2N2 + 1 Grassmann pairs have been used 
up through N2 functions Bn- This means that a single Grassmann pair is missing, hence one 
of the Bn should begin by skipping the Grassmann pair that precedes it. There are N2 ways 
this can happen and the result is that there are now M — 2N2 — 1 Grassmann pairs to be 

distributed among A^^i places. The whole scheme goes on this way. Each time one uses 2N2 
from the first 2N2 + i Grassmann pairs in the N2 Bn functions. The ways this can be done 



are equal to the ways one can distribute i indistinguishable objects at N2 distinguishable 

/ 

places. The result is I " I and has to be multiplied by 




M -2N2 



, which 
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are the ways the remaining Grassmann pairs can be distributed in the A^i places. Summing 
for all i we have the result 

M-NI-2N2 ( ^ + - 1 \ / M -2N2-i\ 

E ■ (17) 

Eq. (17) is true when iV2 7^ 0. When A^2 = the corresponding result is simply the ways we 
can distribute the M Grassmann pairs in A''! places, that is 



(18) 



Since the products coming from the Pexp function are ordered they are all unique, which 
means that each one of them registers directly, i.e. without any weight factor. 

Our next step is to consider how do products of exactly Grassmann pairs at- 
tributed to exp Y^nininCn arisc, if the unused Grassmann pairs from the substructure 
A^i, N2 are 1/3 = M — iVi — 2N2. Obviously, the latter products originate from the component 
of the exponential term 

ITTT n E LAnAn, ■ (19) 

The sum in the above expression is extended over all the Grassmann pairs that have not been 
used in the previous substructure, since the entrance of anyone of the rest M — v-^ pairs in this 
sum will have null contribution to the final result. In short, given that the double appearance 
of a Grassmann variable in any product causes it to vanish, we need to find the number of 
terms that survive. This number is equal to the number of the ways the z/3 Grassmann pairs 
can be distributed in N^, places. Each place, here, corresponds to the position the particular 
Grassmann pair will enter the product. Generally speaking 1/3 Grassmann pairs are allowed 
to occupy the first position. For every Grassmann in the 1st position, 1^3 — 1 Grassmann 
pairs arc allowed to occupy the 2nd, etc. All the allowed products arc therefore the orders of 
z/3 per iV3, that is (z/3) = v-^iyj, — 1) • • • (1/3 — A^3 + 1). But these terms are not all different. 
Two Grassmann pairs, given that they contain two Grassmann variables, may change their 
position without having any effect on their relative sign. So the products of A^3 pairs are 
identical if they contain the same Grassmann variables. We may then collect these identical 
products. Their number coincides with the number of the permutations of N^, objects to As 
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places, that is N^l. The number of different terms from the part of the exponential given by 
(18) diminishes then to combinations of per 

{i^3)ns 1/3(1^3 - 1) •■• (^^3 - + 1) 1/3(1^3 -l)---(z^3-iV3 + l)(l/3-A^3)! 



(20) 



The factor accompanying each of these terms will be then N^l, divided by the same factor 
entering (19) in the denominator. The result will be an overall unit factor. 

Next we ask how do products of exactly 2N4 Grassmann pairs ^nCn enter through 



exp 



J2n^m Cn^nCm^mDnm , if the unuscd Grassmann pairs from the previous substructures 



are 1^4 — M — Ni — 2N2 — N^. These products originate from the part of the exponential 

2 A^4 / 1^4 



j=l \nj=lmjj^nj 



N4I 



(21) 



The number of terms that survive is the number of the ways we can distribute the 1/4 Grass- 
mann pairs to 2N4 positions, since every function Dnm provides two positions of Grassmann 
pairs. The number of allowed products are the orders of 1^4 per 2N4, that is (i^4)2iV4- In 
order to isolate different terms we must divide by the ways the A^4 functions Dnm could 
appear in the product. There are A^4! such ways. But since every Dnm function is equal to 
a Dmn function, we have to divide by a factor of 2 for every such function. We have N4 
such functions, so the total factor by which we divide is 2^"*. Summarising, the number of 
surviving, different terms is 



1 _ (2iV4)! 




(22) 



The factor accompanying each of these terms will then be 2^''A^4!, divided by the A^4! which 
appears in (21), the net result being an overall factor of 2^-*. 

Collecting the above findings wc arrive at the conchision that the total number of terms 
in a structure A^i, N2, N^, N4 that survive, while being different, is 



R{M,Ni,N2) 




M -Ni- 2N2 

N3 



{2N4)\ ( M-N1-2N2-N3 
2^ 2iV4 



, (23) 
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where 



R{M,N,,N2) 




M -2N2-i 



, if A^2 7^ 

■ (24) 



The accompanying factor for every term is 2^*. 

The above considerations are important in the building up of the program. First the 
knowledge of the factor accompanying each term is necessary, since in this way one avoids 
producing identical terms and subsequently summing them up. Also the number of different 
terms existing in each substructure is needed to parameterise the dimensions of the arrays 
that are used in the communication between the subroutines of the program, as it will 
become apparent in the next section. Moreover, the theoretical prediction of the number of 
terms expected in the final result is used to verify the output of the program. This is most 
important especially for large M, where there are numerous surviving terms. 

4. Presenting the structure of the program 

The program has been written in FORTRAN 90. The reason for this choice is the use of 
allocatable arrays in all the subroutines. This enables us to parameterise the dimensions of 
these arrays for different values of M, since the size of some of them grows drastically with 
M. We are also able to manage the available memory effectively. 

According to our previous analysis we shall first produce in the program all the groups 
of integers A^i, N2, N^, N4 that define via eq. (16) the compatible structures for given M. 
This is accomplished in the subroutine STRUCTUREMATRIX, which performs four loops 
at which the integers A^i and A^3 acquire all values from to M and the integers N2 and N4 
all values from to M/2. Once a successful combination for the values of these integers is 
found, according to (16), it is copied to a temporary matrix, for which its column dimension 
is known a priori. This dimension is the number of all the different structures that exist for 
a particular M, a number that grows relatively slowly with M.^ This number can be easily 
evaluated for all M and has been inserted to the data of the program for all M < 20. For 
^For example it equals 486 for M = 20. 
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M > 20 the program calls, where it is needed, subroutine DIMOFSTRUCTUREMATRIX, 
which produces the number of different structures for the given M. Given that we are 
interested in grouping the resulting terms according to their power of T, we group the 
structures according to this power. It is easily checked from eqs. (1), (11) and (13) that this 
power is M — 3 — A^2 — A^4- The final result is inserted into a matrix which is the output of 
STRUCTUREMATRIX. 

The surviving terms in the final result are produced in subroutine MULTIPLYALL. The 
latter first calls subroutine STRUCTUREMATRIX and then for every particular structure it 
calculates first the compatible substructure Ni,N2, followed by the compatible substructure 
A^3 and finally the compatible substructure N4. 

Three complementary routines to MULTIPLYALL have been constructed. The first is 
subroutine COMBINATION ASCEND, which produces all the combinations of the L natural 
numbers from 1 to L in groups of K and places them in ascending order in the lines of a 
matrix which constitutes the output of the subroutine. An array A with K elements is used 
for temporary inputs. The numbers from 1 to X are placed in the first line of the output 
matrix, with the last K-th column being marked as the change column. This means that the 
lines which follow in the output matrix are written with the first K — 1 elements the same 
while the K-th element acquires all the values from + 1 to L. The routine is repeated with 
the change column lowered by one. If the change column is the J-th, it is checked whether 

Aj<M -N + J . (25) 

If (25) holds, the elements with index greater than J acquire the values 

Aj+j = Aj + I, {I = J+1,K). (26) 

Then the array A is placed in the output and K-th column is marked again as the change 
column, repeating the procedure. In the opposite case, where (25) is false, the change column 
is lowered by one, until condition (25) is fulfilled. An example of the output of this subroutine 
is given in Table 1. 

The second subroutine is ORDERS, which finds the orders of L natural numbers from 
1 to L in groups of K and places them in the lines of the output matrix C. The matrix 
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C is filled column by column. In every column two integers determine the filling course: 
The first is the number of identical elements of the previous column. This integer is called 
STEPGREAT and for the i-th column it equals ^^zj^- This expression gives the number 
of the orders of the L — i + 1 numbers which are left after the first i — 1 have correspondingly 
occupied the first i — 1 columns. The second integer, called STEPSMALL, is the number 
of repetitions of the same element in the i-th column and equals STEPGREAT divided by 
L — i + 1. An example of the gradual filling of matrix C is shown in Table 2. 

The third subroutine is COMBINATIONPERTWO which finds the combinations of the 
L' natural numbers from 1 to L' in groups of K' after they have been grouped in couples. 
Each group is placed in ascending order in the output matrix C. Each couple always enters 
as nm with n < m. Matrix C occupies 2K' columns so that each couple is located at adjacent 
columns. To achieve this COMBINATIONPERTWO calls first COMBINATIONASCEND 
(with L = L' and K = K') and stores the result in matrix CA. Then it calls ORDERS 
(with L — L' — K' and K — K') and stores the result in matrix CB. Matrix CA produces 
the odd columns of matrix C, that is the first elements of the couples. For every line of 
CA the remaining L' — K' elements are found and placed in a temporary array B. The 
second element of every couple is found from matrix CB after its elements are read as 
indices of the array B. Because of this last action subroutine ORDERS is only called once in 
COMBINATIONPERTWO. In the production of the group of couples it is checked whether 
the second element is greater than the first. If this is not so the subroutine proceeds to the 
next combination. An example of the output of the subroutine is given in Table 3. 

With the three complementary subroutines that have been described above the calcula- 
tion of the substructures can be carried out. First the calculation of the substructure A^i, A^2 
is produced. By calling GOMBINATIONASGEND from MULTIPLYALL the combinations 
of the positions of the N2 Bn functions from a total of A^i + N2 vacancies are found and 
placed in the matrix GPL. For every such combination, which is represented as a line of 
GPL, one has to produce all the compatible products of functions An and i?„. First wc 
shall assume that the integers n, which characterise the functions An and Bn, run from 1 
to M and that they must appear in an ascending order, from small to large numbers. We 
then allocate two matrices GTEMPl and GTEMP2 with equal dimensions. The number 
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of columns equals A^i + N2 and the number of lines is given by eq. (17) or (18). The last 
value is returned by function C0UNTTERMS2. For every element in the positions from 1 to 
A^i + N2 the maximum and minimum allowed values are found and are placed in the arrays 
IMAX and IMIN respectively. One has to bear in mind that every function consumes 
two consequent integers, so when one such function, characterised by integer n, appears in 
the i-th place, then the integer in the {i + l)-th place cannot be smaller than n + 2. The 
calculation begins by inserting in the first column of matrix CTEMPl and in the different 
lines all the allowed values for the first integer, as they are dictated by the arrays IMAX 
and IMIN. In the second run all the allowed combinations of the first two elements will be 
copied in matrix CTEMP2. To this end, for each hue of CTEMPl the allowed values of 
the next element are determined. These allowed values can range from the maximum to the 
minimum value of the arrays IMAX and IMIN but they must also fulfil additional constraints 
which depend on the preceding integer. If the preceding integer is n, then the next integer 
will be n + 1 when the corresponding function is of the type A and n + 2 if it is of type 
B. So each hne of CTEMPl does not necessarily correspond to the same number of hues 
of CTEMP2. The calculation proceeds with the combinations of the three first elements 
copied back to CTEMPl with the data taken from CTEMP2 and so forth. Once all the 
elements are exhausted, their allowed combinations will have been accommodated in one of 
the matrices CTEMPl or CTEMP2 and are then copied in matrix C12. In the following 
example we suppose that M = 5, A^i = 2, A^2 = l and that the lone function B occupies 
the middle place among the three available for the functions A and B. The array IMIN has, 
subsequently, acquired the values ^1 2 4^ and the array IMAX the values ^235 
The calculation of the allowed combinations proceeds as follows 




CTEMP2 

^ 2^ 

1 3 

2 3 



CTEMPl 

^12 4^ 

12 5 

1 3 5 

2 3 5 



The compatible substructure A^3 for every line in the substructure A^i,A^2 is produced 
with the call of the subroutine COMBINATIONASCEND with L ^ M - Ni - 2N2 and 
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K — Ni and the output is stored in matrix CI. Also, for every line of CPL the remain- 
ing elements which have not been used are determined and placed in the array LEFT34. 
Matrix CI contains the indices of the array LEFT34, i.e. the presence of 2 in a line of CI 
has the meaning of the 2nd element of array LEFT34. With this representation the sub- 
routine COMBINATIONASCEND is called only once from MULTIPLYALL to produce the 
substructure A^3. 

The compatible substructure for every line in the substructures A^i , N2 and N^i is 
produced with the call of the subroutine COMBINATIONPERTWO with V = M - Ni - 
2N2 — N3 and K' — A^4 and the output is stored in matrix C2. Also for every line of CPL and 
every line of CI the remaining elements which have not been used are found and placed in the 
array LEFT4. Matrix C2 contains the indices of the array LEFT4, i.e. the presence of 2 in a 
line of C2 has the meaning of the 2nd element of array LEFT4. With this representation the 
subroutine COMBINATIONPERTWO is called only once from MULTIPLYALL to produce 
the substructure A^4. 

Until now we have described the production of the sequence of functions A, B, C and D 
in the terms that survive in the final result. The indices n which characterise the functions in 
this description are considered as natural numbers from 1 to M and their order is considered 
to be an ascending one. Clearly, the order of these natural numbers is the time order given a 
priori by the user. This time order is stored in the array TOR. So in order to obtain the final 
result we have to read the elements in the structure which has been described as indices of 
the array TOR. The output of MULTIPLYALL is, for every structure, a series of lines which 
represent the products of the functions (10)-(13). Their representation, which is given by 
integers, is listed in Table 4. It should be noted that even though the functions An, and 
Dnm may happen to be represented by the same integer, they cannot be confused, because 
they have a specific order of appearance. Specifically the non-commutative objects of the 
substructure A^i, A^2 appear first, followed by the commutative objects of the structure N2, 
and finally by those of the structure N/^. 

The output of MULTIPLYALL is written in a series of files with the first one designated 
as FOOl. When a file exceeds 18,000,000 terms then the next structure will be written in the 
next file. Each line in a structure is accompanied by a factor which is just 2^*, according to 
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the discussion of the previous section. In the cases for which the substructure Ni — N2 — 
enters, this factor is multiphed by 2. This comes from TriiS^u) — 2 = 4 — 2 = 2. The 
output of MULTIPLYALL is written in the form of integers and so it can be used for further 
calculations by the user. 

The program gives in the output of MULTIPLYALL the products of surviving functions. 
However the Lorentz trace has yet to be performed. This trace involves all the functions 
in the Ni,N2 substructure and this is the reason these objects have been considered as 
being non commutative. The Lorentz indices /i, u in the functions (10), (11) have been 
considered indefinite, that is they have not been defined at this point by the number which 
represents the function. So when groups of these types of functions are found together it is 
assumed that the Lorentz trace has to be performed, according to the exact order by which 
the corresponding functions appear. The following example illustrates the "translation" of 
a typical product of three such functions 

202 301 204 ^ (202)^, (301),^ (204)<,^ ^ (yla)^^ (B,),, (A^),^ . (27) 

In order to perform the Lorentz trace and arrive at the final output we have constructed 
the subroutine TRACE. For the computation of the trace only the part {s^^Pn,^ — £vPn,ti) 
from functions and the part (e^+^e" — e^'^^e^ from functions Bn need to be considered. 
The following representation is used 

el 300 + n, ^ 200 + n , (28) 

where n is a two digit integer. Then the relevant part of the functions A^ and 5„ for the 
trace is represented by two matrices 




300 + n 200 + n 
200 + n 300 + n 



(29) 



eT'eZ - eV\l - " I I I . (30) 



300 + n + 1 300 + n 
300 + n 300 + n + 1 
In a product of functions A^ and where the Lorenz indices are active, the second index 
of one function is saturated with the first index of the next. This continues until the last 
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function in the product is reached, whereupon its second index is saturated by the first 
index of the first function. To accomphsh this in the representation (29)- (30) it suffices to 
perform internal products between the elements of the 2nd column of the 2x2 matrix that 
represents one function and the 1st column of the 2x2 matrix representing the next function. 
The elements of the 2nd column of the last matrix are coupled with the elements of the 1st 
column of the first. The procedure has been appropriately parameterised in the program to 
exhaust all the combinations. 

The rule • p„ = is also taken into account. On the other hand, we have not applied 
the replacement associated with momentum conservation, since it spoils the symmetry of 
the terms belonging to the same structure in addition to increasing the number of terms, 
especially for large M. 

The result of the trace, for each term, corresponds to the sums of products of internal 
products of £'s and p's. These sums arc again represented by series of integers. An internal 
product is represented by a six digit integer, according to the prescription of Table 5. 

The factor of a product is represented by an integer less than or equal to 200000. Such 
factors also signal both the beginning of the product of internal products, as well as, its 
ending. This can be summarised in the next example 

-1 304301 1 304303 302203 301202 204204 -3 304203 303302 ^ 

— ■ 6'^ + ■ ■ ps ■p2 P4-P4 — S ■ ps ■ . (31) 

The powers of 2 and T associated with the functions (10), (11) and (12) are incorporated 
into the factor of every term. The representation of the remaining functions are given through 
integers according to Table 6. 

The output of TRACE is registered in two series of files. The first series begins with file 
GOOl and it contains the result of surviving terms grouped according to the structure they 
belong and arc represented through integers according to Tables 5 and 6. The content of these 
files, as it is in integer form, may be used for further processing by the user. The second series 
begins with file HOOl and it contains the surviving terms in function-like form. The terms 
entering files GOOl and HOOl correspond to the terms stored in files FOOl, and so on. There 
is, also, the option to produce with TRACE only the first term of each structure. In that 
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case TRACE avoids reading through the whole series of terms produced by MULTIPLYALL. 
Instead it calls MULTIPLYALL again and asks for the production of only the first term of 
each structure. MULTIPLYALL writes the result in only one file, FFOOl, which is then read 
by TRACE. With this procedure the timesaving in TRACE is enormous. 

The different structures that exist for the current M are listed in file STRLOG, along 
with information about the files used for the output and the consumed time for the run of 
the program. 

5. Results 

At the beginning of its run the program asks the user to insert the number of external 
gluons M. Because the output of subroutine MULTIPLYALL occupies storage space which 
increases rapidly with M, the user is warned about the approximated space that will be 
required. This space has been measured from the output of the program until some value 
of M. For greater values of M it is extrapolated through a fitting function of the form 
y = exp[a(ln(a;))^ + 6 ln(a;) + c], where y is the storage space, x the total number of surviving 
terms and a, b and c free parameters which have been determined through a fit on the 
existing data. The same function has been used for the extrapolation on the storage space 
required by TRACE. Then the program requests for the specific time order for which the 
calculation will be carried out and a check is performed on the correct input of this time 
order. 

Because the output of subroutine TRACE occupies large storage space the user has the 
option to ask for the full result of this subroutine, or for the calculation of only the first term 
of every structure, or to skip the subroutine entirely. Each choice of output of TRACE is 
accompanied by a warning about the required storage space. 

In Table 7a wc show the result of MULTIPLYALL for M = 2 and in Table 7b the result 
of TRACE in function-like form for the same value of M. The 23 terms for M — 3 are 
shown in function-like form in Table 8. The number of surviving terms grows rapidly when 
M increases (for M — 15 the surviving terms are 430,576,126). 

In Figure 2 we have plotted the number of surviving terms (left axis - curve (a)) against 
M. In the same figure we have plotted the consumed storage space for the output of MULTI- 
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PLYALL (right axis - curve (b)), which is more compact than the full final result of TRACE. 
At M = 15 this space grows to about 22.2 Gbytes, making difficult for us to have the full 
result written on our hard disc for greater M's. 

In Fig. 3 we have plotted the number of structures as a function of M (left axis - curve 
(a)). On the same graph we have depicted the average number of surviving terms per 
structure as a function of M (right axis - curve (b)). It is seen that the number of different 
structures does not grow rapidly, in contrast to the surviving terms per structure. So the 
time needed for the completion of the program is mainly due to the intense growth of number 
of these terms. 

In Fig. 4 we have plotted the consumed time for the completion of the calculation of 
MULTIPLYALL (full line) against M (curve (a)). The consumed time is the most crucial 
parameter since it exhibits the capabilities of our program, especially for large M. This is 
only 2.74 seconds for M = 10 and 1.75 hours for M = 15 (on INTEL CELERON 1.8GHz). 
The capability of the current program is more evident if it is compared with the older program 
we had developed in [13] and where the calculation is carried out in the straightforward way 
without the use of structures. With respect to that program we have made the substitution 
of eq. (12) in order to group functions with the same Grassmann part in the same overall 
function, as the case is for the present program and have recorded the time consumed for the 
calculation without the subroutine TRACE. This allows a comparison to be carried out on 
the same footing. The result is the dotted fine (curve (b)). The slashed line in Fig. 4 (curve 
(c)) is the consumed time in our current program for MULTIPLYALL to be completed but 
with suppressed writing in the output files which occupy enormous storage space. In this 
way we carry calculations for greater values of M. It is also evident, from the comparison 
of curves (a) and (c) that a large fraction of the consumed time for MULTIPLYALL is due 
to the output. In curves (a)-(c) we have used only the time consumed in MULTIPLYALL 
and have not added the time consumed in TRACE because the larger part of the time in 
the latter is consumed in writing the output in the appropriate files and not in the actual 
calculation. The consumed time in TRACE as function of M is presented by curve (i) for 
the full result and by curve (ii) for the result with only one term per structure. 

In Fig. 5 we have plotted the consumed time, as in Fig. 4, but against the surviving 
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terms. The consumed time (full line) in our new program as a function of the surviving 
terms can be fitted through a linear function. The consumed time in our old program grows 
much more rapidly as function of the surviving terms (dotted line). In that case it can be 
fitted by a second grade polynomial which has been extrapolated to larger values of M. This 
fitting procedure leads to a value of about 1910 years for M — 15 and can be compared 
to the 1.75 hours for the same result with our new program. The difference is due to the 
appearance in the intermediate calculations of our old program of a vast number of terms 
which cancel in the final result. 

In the overall result divergent terms enter only for M < 4. These terms have been 
compared with the analytic calculations for the divergent terms for M = 2, M = 3 and 
M — A [10] and have been found in total agreement. Also the number of surviving terms 
according to the program have been compared with the theoretical formulas (23)- (24) and 
have been found equal. The number of surviving terms between the old and the new program 
have also been found in agreement. 

The methods used for the computation of the master formulas do not place an upper 
bound on the value of M for the calculations. In the specific program we have written 
the different representations through integers assuming that M is a two digit integer. The 
limitations are imposed mainly by the storage space and of course, unavoidably, by the 
necessary time which, even for linear growth, will reach high values for increasing M 's. The 
consumed time, though, is expected to shrink proportionally to the increase of the speed of 
the processor to be used. Another limitation arises from the maximum integer that can be 
encoded in the processor one uses for the running of the program. This sets an upper bound 
on the number of the elements that can be allocated to a certain array. The size of the 
matrix that grows more rapidly is the output of COMBINATIONPERTWO which is called 
from MULTIPLYALL to produce the substructure N4. In processors that use 4-byte integers 
overfiow is expected to arise for M > 20, while in processors that use 8-byte integers the 
overflow is expected to appear for M > 34. It is possible to improve these values of M, in 
any case, by defining the respective array in a more complicated way, but we have avoided 
doing so, since the considered values of M are already quite high. 
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6. Concluding remarks 

In this paper a computational algorithm has been presented aiming towards calculation 
of one-gluon loop Feynman diagrams, including ghost loop contributions, with M external 
gluon attachments, on the basis of the master formulas derived in [10]. Compared with 
our previous attempt [13], the reasoning of the calculation has been completely changed, 
resulting to a different algorithm. The basic innovations are the grouping of similar terms 
together and the dramatic reduction of computation time, thus allowing the calculation for 
quite large values of M . All this being said, there remains, of course, the problem of carrying 
out the integrations over the Feynman parameters, a task which is challenging in itself, both 
on the analytical and the numerical front. With respect to the latter, it would be of interest 
to assess the extent to which recent, relevant, considerations [14] can be apphed for the 
succesful completion of the task, at least up to some reasonably high value of M, whose 
Grassmann variable integration part was accomplished in this work. 

Finally, one may express the hope that the particular feature of the constructed algorithm, 
namely the ability to expedite integrations over a multivariable set, a subset of which is 
Grassmannian, could find wider applications to analogous situations that may arise in other 
physical problems wherein Grassmann variables make their entrance. Within the context of 
the present application, it would be of interest to apply the particular algorithm developed in 
this work to the two-gluon loop M — A case, the corresponding master expressions for which 
have been derived in [11]. As a first attempt, one could restrict the relevant computation to 
the divergent term associated with the M — 2 configuration and verify the consistency with 
second order corrections to the running coupling constant in pQCD. 
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Tables 

12 4 
1 3 4 
V 2 3 4 j 

Table 1. The output matrix of subroutine COMBINATIONASCEND with L = A and 

K ^3. 



/ 


1 


0^ 






1 


2 






1 









1 


3 






2 









2 


1 






2 









2 


3 






3 









3 


1 




\ 


3 


J 




I 


3 


2 


/ 



Table 2. The gradual filling of the output matrix of subroutine ORDERS with L — 3 and 

K = 2. 

^13 2 4^ 

14 2 3 
^12 3 4^ 

Table 3. The output matrix of subroutine COMBINATIONPERTWO with M = A and 

Function Representation 
An 200 + n 

Bn 300 + n 

Cn n 

Dnm 100 *n + m 

Table 4- The representation of functions through integers in the output of subroutine 
MULTIPLYALL. The numbers n, m are to be considered as two digit integers. 
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Internal Representation 
Product 



Pn 


■Pm 


200000+n*100+200+m 




■Pm 


300000+n*100+200+m 






300000+n*100+300+m 



Table 5. The representation of internal products through integers in the output of TRACE. 



Function Representation 

300+n 
20000+n*100+m 

30000+n*100+m 
(=30000+m*100+n) 

Table 6. The representation of functions through integers in the output of TRACE. 



M= 2, TIME ORDER: 1 2 






STRUCTURE:( 1) [T**( 


-2)] 




TERMS= 




1. 


FACTOR= 




4. 


102 






STRUCTURE:( 2 0) [T**( 


-1)1 




TERMS= 




1. 


FACTOR= 




2. 


1 2 






STRUCTURE:( 2 0) [T**( 


-1)1 




TERMS= 




1. 


FACTOR= 




1. 


202 201 






TOTAL # OF 






TERMS= 




3. 



Table 7a. The output of MULTIPLYALL for M ^2 and the time order Ui <U2. 



■ e^dndmG{un, Um) 
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M= 2, TIME ORDER: 1 2 



-(pi**2/2)*g**(2 )*d4(pl +p2 )*TrC(tGa2 tGal )* 
inf 1 u2 2 2 

I dT I du2 I dul 8(u2 ,ul ) * Exp{T S S [pn.pmG(un,um)]}* 

n=l m=n+l 

STRUCTURE: ( 1) TERMS= 1. 

2. T**(-2 ) 

{+ el .e2 dl d2 G(ul ,u2 ) } 

STRUCTURE: ( 2 0) TERMS= 1. 

2. T**(-l ) 

{+ S.(i.ne.l )[el .pidl G(ul ,ui)] S.(i.ne.2 )[e2 .pid2 G(u2 ,ui)] } 
STRUCTURE: ( 2 0) TERMS= 1. 

4. T**(-l ) 

{+ 

( +2 c2 .pi el .p2 -2 e2 .el p2 .pi ) } 

TOTAL # of Terms = 3. 

Table 7b. The function like output of TRACE for M = 2 and the time order Ui < U2. The 

Memorandum is similar to Table 8. 
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M= 3, TIME ORDER: 3 12 



-(pi**2/2)*g**(3 )*d4(pl +p2 +p3 )*TrC(tGa3 tGa2 tGal )* 
inf 1 u2 ul 3 3 

I dT I du2 I dul I du3 8(u2 ,ul ,u3 ) * Exp{T S S [pn.pmG(un,um)]}* 

n=l m=n+l 

STRUCTURE:( 11) TERMS= 3. 

2. T**(-l ) 

{+ S.(i.ne.3 )[e3 .pid3 G{u3 ,ui)] el .e2 dl d2 G(ul ,u2 ) 
+ S.(i.nc.l )[cl .pidl G(ul ,ui)] c3 .c2 d3 d2 G(u3 ,u2 ) 
+ S.(i.ne.2 )[e2 .pid2 G(u2 ,ui)] e3 .el d3 dl G(u3 ,ul ) } 
STRUCTURE: ( 110 0) TERMS= 2. 

8. T**(-l ) 

{+ d(ul -u3 ) 

( +2 e3 .e2 el .p2 -2 e3 .p2 e2 .el ) 
+ d(u2 -ul ) 

( +2 e3 .el e2 .p3 -2 e3 .e2 el .p3 ) } 
STRUCTURE: ( 3 0) TERMS= 1. 

2. T**( ) 

{+ S.(i.ne.3 )[e3 .pid3 G(u3 ,ui)] S.(i.ne.l )[el .pidl G(ul ,ui)] S.(i.ne.2 )[e2 .pid2 G(u2 ,ui)] } 
STRUCTURE: ( 2 10) TERMS= 3. 

4. T**( ) 

{+ S.(i.ne.2 )[e2 .pid2 G(u2 ,ui)] 

( +2 e3 .pi el .p3 -2 e3 .el p3 .pi ) 
+ S.(i.ne.l )[el .pidl G(ul ,ui)] 

{ +2 c3 .p2 c2 .p3 -2 e3 .e2 p3 .p2 ) 
+ S.(i.ne.3 )[e3 .pid3 G(u3 ,ui)] 

( +2 e2 .pi el .p2 -2 e2 .el p2 .pi ) } 
STRUCTURE:( 3 0) TERMS= 1. 

-8. T**( ) 

{+ 

( +1 e3 .pi e2 .p3 el .p2 -1 e3 .e2 el .p2 p3 .pi -1 e3 .el e2 .p3 p2 .pi 
+1 e3 .e2 el .p3 p2 .pi -1 e3 .pi e2 .el p3 .p2 +1 e3 .p2 e2 .el p3 .pi 
+1 e3 .el e2 .pi p3 .p2 -1 e3 .p2 e2 .pi el .p3 ) } 



TOTAL # of Terms = 10. 



Memorandum: 



d4(pl+p2+p3 ) 


S^'^^pi +P2 +P3) 


pi**2 




d(u2 -ul ) 


5{U2 — Ul) 


c2.e3 


£2. £3 


TrC(tGa3 tGa2 tCal ) 




inf 


oo 


d2d3G(u2,u3) 


d2d3G{u2,U3) 


cl.p2 


■P2 


8{u2 ,ul ,u3 ) 


e(u2,Ul,U;i) 


S 


E 


T**(-l) 


T-i 


pl.p2 


Pi ■P2 


S.(i.ne.2 )[e2 .pid2 G(u2 ,ui) 


J2i^2^^ ■ Vid2G{u2,Ui) 


1 


/ 











Table 8. The function like output of TRACE for M — 3 and the time order U3 < ui < U2. 
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Figures 







Figure 1 Illustration, for M = 4, of the classes of one-gluon loop Feynman diagram (right 
side of arrow) which are simultaneously accommodated by the corresponding master 

formula depicted on left side of arrow. 
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2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 



M 

Figure 2 (a) The number of surviving terms (left axis) as function of M. (b) The storage 
space occupied by the output of subroutine MULTIPLYALL (right axis). 
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M 

Figure 3 (a) The number of structures as function of M (left axis), (b) The average 
number of surviving terms per structure as function of M (right axis). 
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2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 



M 

Figure 4 Curves (a)-(c) depict the consumed time for the completion of subroutine 
MULTIPLYALL as function of M . (a): Our new program with the output stored in files, 
(b): Our old program [13] with the output stored in files, (c): Our new program with 
suppressed output. Curves (i)-(ii) depict the consumed time for the completion of 
subroutine TRACE as function of M, for our new program, (i): Full result, (ii): Result 
with only the first term of each structure. (Time in all cases in measured using as 
processor INTEL CELERON at 1.8 GHz). 
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10^^ ^1 iiiiiiii I iiiiiiii I iiiiiiii I iiiiiiii I iiiiiiii I iiiiiiii I iiiiiiii I iiiiiiii I iiiiiiii I iiiiiiii I I lira 



10 



10 



10^ 



10^ 



O 10^ V 



C/3 



10^ 



10^ -k- 



10^ ^ 



10^ -= 



10-^ -= 



1 1 




ffi 



10° 10^ 10^ 10^ lO"^ 10^ 10^ lo'^ 10^ 10^ 10^° 10^^ 



# of Surviving Terms 



Figure 5 The consumed time in subroutine MULTIPLYALL ((a)-(c)) and TRACE 
((i)-(ii)), as in Fig. 4, as function of the surviving terms in the final result. Curve (a) has 
been fitted with by a linear function and (b) has been extrapolated up to the surviving 

terms for M = 15 by a quadratic function. 
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